home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / graphic / jpegsrc4.zip / JREVDCT.C < prev    next >
C/C++ Source or Header  |  1992-11-10  |  14KB  |  370 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on an algorithm described in
  11.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  12.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  13.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  14.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  15.  * We use their alternate method with 12 multiplies and 32 adds.
  16.  * The advantage of this method is that no data path contains more than one
  17.  * multiplication; this allows a very simple and accurate implementation in
  18.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  19.  */
  20.  
  21. #include "jinclude.h"
  22.  
  23. /*
  24.  * This routine is specialized to the case DCTSIZE = 8.
  25.  */
  26.  
  27. #if DCTSIZE != 8
  28.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  29. #endif
  30.  
  31.  
  32. /*
  33.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  34.  * on each column.  Direct algorithms are also available, but they are
  35.  * much more complex and seem not to be any faster when reduced to code.
  36.  *
  37.  * The poop on this scaling stuff is as follows:
  38.  *
  39.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  40.  * larger than the true IDCT outputs.  The final outputs are therefore
  41.  * a factor of N larger than desired; since N=8 this can be cured by
  42.  * a simple right shift at the end of the algorithm.  The advantage of
  43.  * this arrangement is that we save two multiplications per 1-D IDCT,
  44.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  45.  *
  46.  * We have to do addition and subtraction of the integer inputs, which
  47.  * is no problem, and multiplication by fractional constants, which is
  48.  * a problem to do in integer arithmetic.  We multiply all the constants
  49.  * by CONST_SCALE and convert them to integer constants (thus retaining
  50.  * CONST_BITS bits of precision in the constants).  After doing a
  51.  * multiplication we have to divide the product by CONST_SCALE, with proper
  52.  * rounding, to produce the correct output.  This division can be done
  53.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  54.  * as long as possible so that partial sums can be added together with
  55.  * full fractional precision.
  56.  *
  57.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  58.  * they are represented to better-than-integral precision.  These outputs
  59.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  60.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  61.  * intermediate INT32 array would be needed.)
  62.  *
  63.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  64.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  65.  * shows that the values given below are the most effective.
  66.  */
  67.  
  68. #ifdef EIGHT_BIT_SAMPLES
  69. #define CONST_BITS  13
  70. #define PASS1_BITS  2
  71. #else
  72. #define CONST_BITS  13
  73. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  74. #endif
  75.  
  76. #define ONE    ((INT32) 1)
  77.  
  78. #define CONST_SCALE (ONE << CONST_BITS)
  79.  
  80. /* Convert a positive real constant to an integer scaled by CONST_SCALE. */
  81.  
  82. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  83.  
  84. /* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
  85.  * causing a lot of useless floating-point operations at run time.
  86.  * To get around this we use the following pre-calculated constants.
  87.  * If you change CONST_BITS you may want to add appropriate values.
  88.  * (With a reasonable C compiler, you can just rely on the FIX() macro...)
  89.  */
  90.  
  91. #if CONST_BITS == 13
  92. #define FIX_0_298631336  ((INT32)  2446)    /* FIX(0.298631336) */
  93. #define FIX_0_390180644  ((INT32)  3196)    /* FIX(0.390180644) */
  94. #define FIX_0_541196100  ((INT32)  4433)    /* FIX(0.541196100) */
  95. #define FIX_0_765366865  ((INT32)  6270)    /* FIX(0.765366865) */
  96. #define FIX_0_899976223  ((INT32)  7373)    /* FIX(0.899976223) */
  97. #define FIX_1_175875602  ((INT32)  9633)    /* FIX(1.175875602) */
  98. #define FIX_1_501321110  ((INT32)  12299)    /* FIX(1.501321110) */
  99. #define FIX_1_847759065  ((INT32)  15137)    /* FIX(1.847759065) */
  100. #define FIX_1_961570560  ((INT32)  16069)    /* FIX(1.961570560) */
  101. #define FIX_2_053119869  ((INT32)  16819)    /* FIX(2.053119869) */
  102. #define FIX_2_562915447  ((INT32)  20995)    /* FIX(2.562915447) */
  103. #define FIX_3_072711026  ((INT32)  25172)    /* FIX(3.072711026) */
  104. #else
  105. #define FIX_0_298631336  FIX(0.298631336)
  106. #define FIX_0_390180644  FIX(0.390180644)
  107. #define FIX_0_541196100  FIX(0.541196100)
  108. #define FIX_0_765366865  FIX(0.765366865)
  109. #define FIX_0_899976223  FIX(0.899976223)
  110. #define FIX_1_175875602  FIX(1.175875602)
  111. #define FIX_1_501321110  FIX(1.501321110)
  112. #define FIX_1_847759065  FIX(1.847759065)
  113. #define FIX_1_961570560  FIX(1.961570560)
  114. #define FIX_2_053119869  FIX(2.053119869)
  115. #define FIX_2_562915447  FIX(2.562915447)
  116. #define FIX_3_072711026  FIX(3.072711026)
  117. #endif
  118.  
  119.  
  120. /* Descale and correctly round an INT32 value that's scaled by N bits.
  121.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  122.  * the fudge factor is correct for either sign of X.
  123.  */
  124.  
  125. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  126.  
  127. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  128.  * For 8-bit samples with the recommended scaling, all the variable
  129.  * and constant values involved are no more than 16 bits wide, so a
  130.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  131.  * this provides a useful speedup on many machines.
  132.  * There is no way to specify a 16x16->32 multiply in portable C, but
  133.  * some C compilers will do the right thing if you provide the correct
  134.  * combination of casts.
  135.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  136.  */
  137.  
  138. #ifdef EIGHT_BIT_SAMPLES
  139. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  140. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  141. #endif
  142. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  143. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  144. #endif
  145. #endif
  146.  
  147. #ifndef MULTIPLY        /* default definition */
  148. #define MULTIPLY(var,const)  ((var) * (const))
  149. #endif
  150.  
  151.  
  152. /*
  153.  * Perform the inverse DCT on one block of coefficients.
  154.  */
  155.  
  156. GLOBAL void
  157. j_rev_dct (DCTBLOCK data)
  158. {
  159.   INT32 tmp0, tmp1, tmp2, tmp3;
  160.   INT32 tmp10, tmp11, tmp12, tmp13;
  161.   INT32 z1, z2, z3, z4, z5;
  162.   register DCTELEM *dataptr;
  163.   int rowctr;
  164.   SHIFT_TEMPS
  165.  
  166.   /* Pass 1: process rows. */
  167.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  168.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  169.  
  170.   dataptr = data;
  171.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  172.     /* Due to quantization, we will usually find that many of the input
  173.      * coefficients are zero, especially the AC terms.  We can exploit this
  174.      * by short-circuiting the IDCT calculation for any row in which all
  175.      * the AC terms are zero.  In that case each output is equal to the
  176.      * DC coefficient (with scale factor as needed).
  177.      * With typical images and quantization tables, half or more of the
  178.      * row DCT calculations can be simplified this way.
  179.      */
  180.  
  181.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  182.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  183.       /* AC terms all zero */
  184.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  185.       
  186.       dataptr[0] = dcval;
  187.       dataptr[1] = dcval;
  188.       dataptr[2] = dcval;
  189.       dataptr[3] = dcval;
  190.       dataptr[4] = dcval;
  191.       dataptr[5] = dcval;
  192.       dataptr[6] = dcval;
  193.       dataptr[7] = dcval;
  194.       
  195.       dataptr += DCTSIZE;    /* advance pointer to next row */
  196.       continue;
  197.     }
  198.  
  199.     /* Even part: reverse the even part of the forward DCT. */
  200.     /* The rotator is sqrt(2)*c(-6). */
  201.  
  202.     z2 = (INT32) dataptr[2];
  203.     z3 = (INT32) dataptr[6];
  204.  
  205.     z1 = MULTIPLY(z2 + z3, F